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The differences are discussed between various next-to-leading order prescriptions for the QCD evolution of parton 
densities and structure functions. Their quantitative impact is understood to an accuracy of 0.02%. The uncer- 
tainties due to the freedom to choose the renormalization and factorization scales are studied. The quantitative 
consequences of the different uncertainties on the extraction of the strong coupling constant from scaling 
violations in deep-inelastic scattering are estimated for the kinematic regime accessible at HERA. 



1. Introduction 

Deeply inelastic scattering (DIS) has a long and 
successful history in the discovery and analysis 
of the substructure of the nucleon. It provides 
one of the cleanest tests of QCD. With the ad- 
vent of HERA, the kinematic range has been ex- 
tended to values up to about lO^GeV^. More- 
over in the Bj or ken- variable x, the range down 
to values of a; ~ 10~^ is now probed in the deep- 
inelastic regime. This extended kinematic cover- 
age allows for detailed measurements of the scal- 
ing violations of the structure function F2{x, Q^) 
and, consequently, of the strong coupling con- 
stant Us- The statistical and systematic exper- 
imental errors for such analyses have been esti- 
mated in refs. 

The determination of as in analyses of DIS 
data from high statistics experiments relies on 
perturbative QCD expansions, implemented in 
complicated fitting procedures. Presently the 
necessary theoretical ingredients are fully known 
only up to next-to- leading order (NLO). In this 
paper, we analyze the theoretical ambiguities 
arising at this level of accuracy, and study their 
implications for the extraction of q;<;(M|) from 
the scaling violation of the structure function F2 ■ 

*Talk presented by S. Riemersma 



The paper is organized as follows. In Section 2 
we give a detailed discussion of the various ap- 
proximations used in different solutions of the 
evolution equations. We examine their quantita- 
tive impact on the predicted scaling violations of 
the parton densities and structure functions. Sec- 
tion 3 is devoted to the uncertainties originating 
from the choice of the renormalization and mass 
factorization scales at this order of the pertur- 
bative expansion. The implications of all these 
effects for the theoretical error Aq;s(A/|) from 
QCD analyses of DIS scaling violations in the 
kinematic range at HERA are presented in Sec- 
tion 4. Section 5 contains the conclusions. 

2. Structure function evolution in NLO 

In this section the details of various representa- 
tions of the evolution equations in NLO and the 
renormalization group equation determining as 
are compared. To keep the notation as compact 
and transparent as possible, the evolution equa- 
tions will be written in a generic manner, which 
covers the non-singlet equations as well as the 
coupled singlet quark and gluon evolution equa- 
tions. Explicit solutions written down in Sections 
2.2 and 2.3 will apply literally only to the non- 
singlet cases. A numerical comparison of the dif- 
ferent prescriptions is performed in Section 2.4. 
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2.1. Evolution equations for parton densi- 
ties and as{Q'^) 

The evolution equations for the (twist-2) parton 
distributions f{x, M^) of the nucleon are given by 

df{x,M^) 



91nM2 



[as{M^) Po{x) + aliM^) Pi{x) 
+ Oial)]^fix,M^). (1) 



As usual X stands for the nucleon's momentum 
fraction carried by the partons, and (g) denotes 
the Mellin convolution, 



A{x) (g) B{x) 



(2) 



M represents the mass factorization scale. In 
eq. (|l|) we have already inserted the perturbative 
expansion of the splitting functions P{x,M^) in 
powers of the strong coupling constant as{M'^) = 
47ras(Af2). Only the first two expansion coeffi- 
cients Po{x) and Pi{x) are completely known so 
far. Hence the solution of the evolution equations 
is presently possible up to NLO. 

To this accuracy the scale dependence of 
as{R^) is governed by 



= ^(3oal{R')~PialiR^) + Oiat), (3) 



aini?2 

with the first two (scheme independent) coeffi- 
cients of the QCD beta function given by /3o = 
II - 2Nf/5 and /3i = 102 - 38iV//3. Here Nf de- 
notes the number of quark flavours. As already 
indicated in eq. (|l|), the renormalization scale R 
is identified with the factorization scale in this 
section, R = M. For the case that R and M are 
chosen to be unequal see Section 3. 

Like the solution of the parton evolution equa- 
tion discussed in detail below, also the NLO 
running of Os can be treated in several ways, dif- 
fering beyond the accuracy of the present approx- 
imation. First, eq. (^) can be simply solved by 
numerical iteration starting from a reference scale 
Rq, or equivalently by using its integrated implicit 
form 
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1 



c^siRl) 



/3oln(|^ 



(4) 



(h f a,iR^)[(3o + PiasiRl)] 
/3o ""{a.iRDiPo + PiasiR^)] 



Second, since in lowest order l/us ^ \n{R'^/Rl), 
one can treat eq. d) in the sense of a power ex- 
pansion in l/\n{R^/RQ), and keep only the pow- 
ers to which the terms in eq. (|^) beyond /?i do 
not contribute. Introducing the QCD scale pa- 
rameter A instead of the reference value as(i?o), 
the result can be written as 



/3i In [ln(i?VA2)] 



/3o ln(i?2/A2) 0i In2(i?2/A2) 
0(ln-''^(i?VA2)) . (5) 



Both these treatments, as well as slightly different 
analytic approximations e.g. by truncating in 
I/fls instead of ag as in eq. have been widely 
used in the context of parton evolution. 

Let us now return to the evolution equation (0). 
This system of coupled integro-differential equa- 
tions has been treated in various manners. Aside 
from expansions in orthogonal polynomials 
which will not be discussed here, the principal 
options are to solve these equations in x-space or 
to consider their transformation using Mellin- A'^ 
moments, 



A'' = / dxx'^-'A{x) 
Jo 



(6) 



In the latter case the convolution (|^) reduces to 
the product 



[A{x)^B{x)f ^A^B^, 



(7) 



and eq. becomes a system of ordinary differ- 
ential equations at fixed N. 

In many analyses, eq. (|^) has been solved nu- 
merically in x-space, see e.g. [D-0, without fur- 
ther reference to the power-series structure in Os . 
We will denote this approach, which is a direct so- 
lution of the renormalization group equations in 
the case of mass factorization, as prescription (A) 
in the numerical comparisons in Section 2.4. On 
the other hand, the transformation of the evolu- 
tion equations to A^-space allows for further 
analytic developments to which we now turn. 

2.2. Analytic solutions and approximations 

The first step towards an analytic solution of the 
evolution equations is to rewrite eq. ([^) in terms 
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of a. 



9a, 



rK)-(8) 



To simplify the notation, the argument of 
is suppressed throughout this section. This for- 
mulation is still fully equivalent to eq. in com- 
bination with eq. (Q). The analytic solution of 
eq. (|^), which is possible in a closed form for the 
non-singlet cases, or a numerical iterative treat- 
ment, can be used to cross-check the numerical 
accuracy of concrete x- and A^-space implemen- 
tations, see ref. jiol] . 

Expanding the r.h.s. of eq. into a power 
series in a,, one arrives at 



or {as) 
das 
1 



(9) 



Poas Po 



The derivation of eq. (|^) consists of two steps: 
first one assumes that the beta function is trun- 
cated for l/tts [and not for a, as in eq. (^] , which 
is of course equivalent to the order considered. 
Then the term arising in the square bracket is 
discarded, since at the same order in as also the 
presently unknown next-to-next-to-leading order 
(NNLO) splitting functions enter. Hence this 
term is beyond the order being considered here. 
The direct solution of eq. will be referred to 
as prescription (B) pj] ] in the numerical compar- 
isons below. 

Keeping the power-series character of the evo- 
lution equations, the final step of the solution 
then leads to 



1 - 



"° ^ /~(ao) 
as ' 



+ 0{al) 
(10) 



in the non-singlet cases, with ao = asiM^). 
Here Mq is the reference scale at which the 
non-perturbative initial distributions are speci- 
fied. With respect to the as expansion, the no- 
tationally more involved singlet-matrix solution 
behaves as eq. (p^, sec ref. As compared 

to this final result, an iterative solution of eq. (|^) 



generates also higher powers of the , the con- 
tribution of which is again beyond NLO in the 
sense of an expansion in powers of a, . In Section 



2.4, we will refer to the solution (10), which has 
practically been employed, for example, in refs. 
H,^, as prescription (C). 

The Mellin-moment technique is a very pow- 
erful tool, allowing the implementation of various 
kinds of approximations. As we will mainly rely 
upon this method in our subsequent numerical 
work, a short reminder about the transformation 
of the results back to x-space is in order. The 
inverse Mellin transformation is performed by a 
contour integral in the complex TV-plane, 

f{x,as)=- dzlm[e*^a;-c^/^=^(a,)], (11) 

7^ Jo 

where C — c + ze*"^. Here all functions en- 
countered satisfy (/^)* — ■ The parameter c 
is chosen about one to two units to the right of 
the rightmost singularity of - at = 1 (0) 
for singlet (non-singlet) cases. All singularities 
are situated on the real axis for the NLO evolu- 
tion, and (f) ~ 37r/4 can be safely used even 
down to extremely low x. The choice oi (j) > 7r/2 
leads to a faster convergence of the integral ( [ll| ) 
as z ^ oo. Using a chain of Gauss quadratures, 
a numerical accuracy better than 10~^ is easily 
achieved. 

2.3. Structure function evolution near Qq 

The structure function F2 is obtained, for NP = 
E? = Q^, from the (anti-) quark and gluon densi- 
ties fq and fg considered in the previous sections 

by 

P2(x,Q')= ^ c,[x,a,(g2)]®/^(a,,g2)_ (^3) 

Here Cq^g[x , as{Q'^)] denote the NLO coefficient 
functions for quarks and gluons 

Ci[a;,as((3^)] = CQs5iq5(l - x) + a,((5^) ci,i(x) , 

(13) 

where we identify the factorization scale M by 
setting A'P ~ . If is not too far from the 
reference scale Qq, i.e. if the evolution distance 
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L = ln(Q^/(3Q) is small, one may Taylor expand 
eq. ( p^ around Qq. In this way ^^2(2;, Q^) can be 
expressed in terms of L, the expansion coefficients 
of the coefficient and splitting functions, Ci(x) and 
Pkiix), the parton densities fi{x,Ql) = fuo{x), 
and the strong coupling as((5§) = oq at the initial 
scale. 

The explicit result for this case will also be 
given for the notationally simpler non-singlet 
structure functions only. We will again switch to 
Mellin moments below. The approximation un- 
der consideration is most easily derived directly 
from eqs. (|l|) and (Q) without reference to their 
solutions discussed above. First the Cs-equation 
(^) is solved up to orders Oq and L^, resulting in 



(14) 



The expansion of the r.h.s. of eq. (12) in connec 
tion with eq. 
0{al) and 0{L 



is then, including terms up to 
, given by 



(15) 



N 



pN N 
'^l,q 



N^2 



/3oPo^)^^]}/< 



One may express this solution in terms of the co- 
efficient function to order a^, Cq{ao,Q^ /Qq) ||1^ 
as 



F^iQ') ^ {(l + aocf:,)+ [cf(ao,QVQo) 

-Cf(ao,l)]}/o^. (16) 

Here the general form of the coefficient function 

Cf(a„gV^^') = co,, + a.Ci^,(QVM2) 

+ alC^^^iQyM^) + Oial) 



(17) 



was used choosing = Qq. The local repre- 
sentation ( p^ ) of the structure function evolution 
will be denoted as prescription (D) in the 
discussion below. 

2.4. Numerical comparisons 

We now turn to the comparisons of the numeri- 
cal results obtained by the prescriptions (A)— (D) 



defined in the previous sections. All these ver- 
sions are equivalent to NLO. Their comparison 
gives information about the uncertainty induced 
at the level of possible approximations. The un- 
certainties due to different choices of the renor- 
malization and factorization scales are discussed 
in Section 3. 

The results will be shown for initial distribu- 
tions which, although representing a somewhat 
simplified input, incorporate all features relevant 
to this study in a sufficiently realistic way. Specif- 
ically, we take in the MS factorization scheme at 
= 4 GeV^: 

xuy = AuX°-^{l - x)^, xdy = Adx"-^{1 ~ x)'^, 
xS = E — xUy — xdy = Asx^^''^{l — xY , 
xg = AgX^°-'^{l - x)^ , xc = xc = . 

(18) 

The evolution is performed for four massless 
flavours, using A^(A^/ = 4) = 250 MeV in eq. 
(^). The SU(3)-symmetric sea is assumed to 
carry 15% of the nucleon's momentum at the in- 
put scale, and the remaining coefficients are fixed 
by sum rules. Also the results on F2 will be shown 
employing — B? = . The massless charm 
evolution does not yield a correct representation 
of F2^'^(a:, Q^). However, the conclusions of this in- 
vestigation are not substantially affected by this 
simplification. 

In Figure 1 the singlet quark and gluon den- 
sities obtained from the solutions (A)— (C) are 
compared, after an evolution to = 100 GeV^ 
and 10^ GeV^. To display the differences clearly, 
the curves have been normalized to the values ob- 
tained for case (C). Between the input scale and 
the highest scale considered in the figure, xE and 
xg increase (decrease) by more than an order of 
magnitude at very low (high) values of x, respec- 
tively. 

The difference between the results based on 



eq. 



(B) and (C), is very small for 10 



< 



X < 0.9, and reaches at most the order of 1% for 
smaller x. This small difference is due to the it- 
erated Pi terms in the solution of eq. (^ . On the 
other hand, the deviations of the evolution (A) 
based upon eq. (|^) from these results are rather 
large. These offsets are not related to any numer- 



5 



1.02 
1 

0.98 
0.96 
0.94 h' 



-1 1 I I I I 1 11 1 1 I I I I 1 11 1 1 I I I I III 1 1 I I I I 1 1 



XL® / xL^^^ 



_i I I 



_i I I I I I I 



1.04 - 



1.02 - 



0.98 



0.96 - 



-| 1 I I I I 1 11 1 1 I I I I 1 11 1 1 I I I I III 1 1 I I I I 1 1 



xg«/xg(^^ 



_I ' 



(A) , M^= 100(10'*)GeV^ 

(B) , M^= 100(10'*)GeV' 
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Figure 1: A comparison of the singlet quark and 
gluon densities, and xg, as evolved using the pre- 
scriptions (A, B) discussed in the text, normalized 
to (C) for the initial conditions ([T^). The results 
are displayed at two representative choices of the 
factorization and renormalization scale M . 
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10 
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Figure 2: A comparison of the evolutions of 
the proton structure function F2(x, Q^) as obtained 
using the prescriptions (A, B, D), normalized to 
(C) close to the reference scale Qq = 4 GeV^. 



leal inaccuracies, since tliose are under control to 



the level of 0.02% |10|. Instead these differences 
are due to both of the two expansion steps dis- 
cussed below eq. with the first one, the re- 
expansion of the beta function, being somewhat 
more important numerically. As will be shown 
in Section 4, these differences in F2 result in a 
shift of as(M|) by 0.003 under the conditions at 
HERA. 

Figure 2 shows how the offset between the so- 
lutions (A) and (B,C) develops, for three typical 
values of x, with increasing in the structure 



function F2. The result of the local represen- 
tation (D), representing an expansion valid for 
L — \n{Q^ /Qq) ^ 1, is also given here. For small 
values of L it agrees with scheme (A). The figure 
indicates that the difference, observed in the pre- 
vious comparison, builds up very quickly in the 
vicinity of Qq, and then stays at about the same 
level above « 20 GeV^. This behaviour is 
readily explained by the large size of as at scales 
around 4 GeV^ : changes by almost a factor of 
two between 4 and 20 GeV^. 
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3. Renormalization and Factorization Scale 
Uncertainties 

The renormalization scale dependence of non- 
singlet parton densities is determined by 



ainM2 



+as{R')f3oPi' In ( ^ 



i?2 



-a,(i?2X (19) 



Here the values of at the factorization scale 
and renormalization scale are, to NLO, 
related by 



as{M^) ^ as{R^) 



l+/3oas(i?')ln 



fR^ 



.(20) 



In fact, one can resum the terms in front of the 
LO splitting functions, , in eq. ( p^ using the 
relation (po|). Thus to NLO accuracy the parton 
densities can be expressed as a function of a single 
scale. 

The structure functions [Q"^) are repre- 
sented by the solution of eq. (|l^) convoluted with 
the coefficient functions according to 



F^(Q2)=[l + a,(i?2)c?^^]/^^(g^i?2) 



(21) 



For the study of the renormalization scale de- 
pendence we identify the mass factorization scale 
M2 = Q2 

The factorization scale dependence of the struc- 
ture function is described by 



i^^(Q^) 



1 



(22) 

)]r(M2), 



where is the solution of 
91nAf2 



(23) 



Here the renormalization scale is fixed by R^ = 
M2. 

Eqs. (|2l|) and (|2|) are independent of the 
choices of the scales R and M, respectively, up 
to NLO. A further improvement of the scale sta- 
bility, up to 0{a^g), can be obtained extending 



the above analysis to 3~loop order, if the corre- 
sponding splitting functions become available. To 
derive an estimate of the quantitative importance 
of the present scale dependences, R^ and will 
be varied in the range between and 4Q2 jn 

Section 4. 




Figure 3: The kinematic range for DIS measure- 
ments at HERA used for the determination of the 
theoretical uncertainties of as (Af|) in Section 4. 



4. Estimation of the Theoretical Error 

Aaa(M|) 

The necessity of using perturbative solutions of 
the renormalization group equations in represent- 
ing the observables implies theoretical errors in 
the determination of as in all experimental anal- 
yses. It is convenient to compare them at a com- 
mon reference scale for which we choose /i^ = 
A/f . Furthermore we refer to the MS renormal- 
ization scheme below. Recent as measurements 
were compiled in ||l^,|l^ . The most precise results 
obtained show that the values of as{M'^) as de- 
termined in deep inelastic scattering and in e^e~ 
annihilation seem to differ. The central value de- 
termined by the DIS data is |l6| 



as{Ml) ^ 0.112 ±0.004, 



whereas in the e^e experiments 
as{Ml) = 0.121 ±0.004 



(24) 



(25) 
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is obtained. Future high statistics measurements 
of as from deep inelastic scattering at HERA may 
help to resolve this difference. For the interpre- 
tation of the data, a careful treatment of the the- 
oretical errors is required. In the following we 
study possible sources of theoretical uncertainty 
in the framework of the NLO evolution. They 
include: 

1) The effect arising from the different repre- 
sentations of Us given in eqs. (^) and (||). 

2) The offsets originating in the different NLO 
prescriptions (A)— (C), respectively based 
upon eqs. ([l]), ^ and (]To|), for the structure 
function evolution. 

3) The theoretical uncertainties due to the 
freedom of choice for the renormalization 
and mass factorization scales. 

To gauge the effect of each source of uncertainty, 
a reference data set was constructed for which 
the different displacements were measured using 
the method in the kinematic regime of HERA. 
The different cuts used are illustrated in figure 3. 
The value of as(Af^) was taken to be 0.112 fixing 
Nf = 4 in the entire range of Q^, taking prescrip- 
tion (A) for the evolution, and using eq. (||) for 
the description of as{Q'^). 

As outlined in Section 2 the solution of the evo- 
lution equations either using eq. (|l]) or eq. (|^) 
lead to different results. In the determination of 
as{M^) a relative shift about 0.003 is implied 
comparing both approaches. A similar differ- 
ence has been observed between various fits in 
the analysis of the BCDMS F^^ data 0, see 
Table 1. 

The use of either the representation for as 
given in eq. or the complete solution induces a 
shift of Aas{M'^) of 0.001 or less. A much smaller 
difference of 0.0001 in as{M'^) is obtained com- 
paring fits using a fully iterative solution (^) of 
the evolution equation with one based on the ex- 
pansion (p^). 

The largest contributions to the theoretical er- 
ror of as(M§) are due to the renormalization and 
factorization scale uncertainties. In the present 
analysis, we vary each scale separately in the 



range to AQ^. The effects are particularly 

strong if no further cuts on are applied, see 
Table 2. If only the range above 50 GeV^ 
is considered, the factorization scale variation in 
the above range results in Aas(M|)|M = ±0.003, 
while the corresponding error for the renormaliza- 
tion scale variation is estimated to be 
r2\\ _ +0.004 
-0.006 ■ 



Aa,(M|)|fl 



5. Conclusions 

Various prescriptions for the NLO evolution of 
DIS structure functions have been discussed, and 
their quantitative differences were studied. The 
resulting i/ieoretica/ uncertainties are rather large. 
At cc ~ 10""*, e.g., the spread reaches about 6% 
for the proton structure function F2{x,Q^) for 
identical inputs. A shift of the central value of 
as(M|) by about 0.003 is observed, when the 
same data in the HERA kinematic region are fit- 
ted with programs using the most differing pre- 
scriptions. This is to be contrasted to the off- 
sets induced by the various TiMTTierica^ procedures. 
They are under control between five independent 
implementations [|[0,|,jlli|l|| of F2ix,Q^) at the 
level of 0.02%, sec ref. 

The dependence of the evolution of F2{x,Q^) 
on the choice of the renormalization and mass 
factorization scales R and M has been inves- 
tigated. By varying both scales independently 
within the range (9^/4 to 4(5^, the theoretical er- 
ror Aas(Af|) of future high statistics Us mea- 
surements from F2 scaling violations at HERA 
has been estimated. If a cut of 50 GeV^ is 
applied, at NLO it amounts to 



Aa,(M|) 



+0.004 
-0.006 



-0.003 
+0.003 



M 



A significant reduction of the theoretical uncer- 
tainties can be expected if the presently unknown 
3-loop splitting functions become available. 
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Program 


AMs/MeV, Nj = 4 


as{Ml) 


xVd.o.f. 


Krivokhizhin et al. Jl8] 


231 ± 20 


0.1115 ±0.0015 


176/152 


Virchaux, Ouraou H 


235 ± 20 


0.1118 ±0.0015 


178/150 


Gonzalez- Arroyo et al. 


220 ± 20 


0.1107 ±0.0016 


180/151 


Abbott et al. [|ol 


233 ± 20 


0.1117 ±0.0016 


198/151 


Furmanski, Pctronzio [pT|] 


270 ± 25 


0.1144 ±0.0016 


181/152 



Table 1: Comparison of the QCD scale parameter Aj^ obtained from NLO fits to the BCDIVIS F^^ data 
employing different evolution programs, taken from ref. ||l^. The resulting values for as(Aff) have been 
computed, using eq. ^ and a bottom threshold of TOf, = 4.5 GeV. 



Cut (GeV2) 


AP = QV4 


M2 = Q72 


M2 = 2Q2 


M2 = 4Q2 


4 
20 
50 


+0.0067 
+0.0032 


+0.0029 
+0.0015 


-0.0068 
-0.0024 
-0.0015 


-0.012 
-0.0044 
-0.0029 


Cut (GeV2) 




= Q2/2 


i?2 ^ 2Q2 


i?2 ^ 4g2 


4 
20 

50 


-0.0076 
-0.0067 
-0.0061 


-0.0037 
-0.0032 
-0.0028 


+0.0032 
+0.0027 
+0.0023 


+0.0059 
+0.0049 
+0.0042 


Cut (GeV2) 


(B), eq. (I) 


(C), eq. ® 


as, eq. (|) 




4 
20 
50 


-0.0027 
-0.0025 
-0.0023 


-0.0001 
-0.0001 
-0.0001 


+0.0010 
+0.0009 
+0.0008 





Table 2: The theoretical shifts on q;s(A/|), from scale variations as well as from different NLO evolution 
and as prescriptions, for F2 data in the HERA kinematic range (see Figure 3). The reference data set was 
generated using prescription (A) for the evolution equations, eq. (^) for as, and imposing NP = = . 
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